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Abstract 

We present an approach to calculate ballistic phonon transport that combines the atomistic 
Green's function (AGF) method with ab initio results. For the inter atomic potential we use 
the harmonic approach. The equilibrium positions of the atoms and the inter atomic force con- 
stants(ifcs) are calculated using the ABINIT program package jl|], which is based on density func- 
tional theory. Therefore, the presented approach is parameter free. From the Green's function 
of the system we determine the density of states as well as the transmission function. The ther- 
mal conductance is obtained within the linear response regime. We apply this approach to bulk 
ZnO and bulk ZnS. Transmission functions for different transport direction for each material are 
presented. A comparison of the transmission function shows, that a ZnO/ZnS interface could be 
a promising phonon blocker. Adding such interfaces in ZnO or ZnS based thermoelectric devices 
could therefore increase the figure of merit. 

PACS numbers: 72.20.Pa, 73.40.-c, 71.55.Gs 
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I. INTRODUCTION 



For thermoelectric devices the figure of merit ZT is the most important characteristic 
number since it is directly connected to the theoretical maximum achievable efficiency of 
the device. The figure of merit is defined as 

ZT=^^ (1) 

where a is the electric conductivity, S is the Seebeck coefficient, T is the temperature, 
is the electric part of the thermal conductivity and Kph is the phonon part of the thermal 
conductivity. In 2001 Venkatasubramanian et a/.Q showed that the figure of merit in ther- 
moelectric devices which are based on super lattice structures can be rather large. The large 
figure of merit is explained by a strong decrease in the phonon heat conductivity due to 
additional scattering of phonons at the interfaces and simultaneously an almost uneffected 
electronic transport across these interfaces. Venkatasubramanian et al. investigated super 
lattice structures that are based on Bi2Te3/ Sb2Te3. Recent experimental results show that 
also ZnO/ZnS based super lattice structures and micro structures are a promising candidate 
as thermoelectric materials js], Q| . It can be grown as sputtered thin films as well as fully 
epitaxially using molecular beam epitaxy. The material parameters can be tuned by doping. 
In this paper we introduce a scheme that combines the fiexibility of the atomistic Green's 
function (AGF) method and the ab initio character of DFT to calculate ballistic phonon 
transport. We use this scheme to investigate the thermal conductance of pure ZnO and pure 
ZnS. ZnO has wurtzite structure while ZnS can occur in wurtzite structure and zincblende 
structure. We show that ZnO/ZnS interfaces are promising phonon blockers regardless of 
the ZnS structure. ZnO has a band gap of 3.37eV and shows an intrinsic n-doping, which is 
related to intrinsic Zn interstitials jsj] . The space group of wurtzite structure is PG^mc. Since 
ZnO in wurtzite structure has four atoms in the unit cell the phonon- dispersion relation has 
12 branches. At the F point the 12 branches divide up into 2A1+2B1+2E1+2E2, where 
both El and both E2 modes are double degenerated. In principle, the phonon dispersion 
relation can be measured by means of inelastic neutron-scattering (ISN) or for the Raman 
active modes with Raman scattering. For ZnO several experimental data with both meth- 
ods have been reportedj6-9|. ZnS has a slightly higher band gap than ZnO of 3.8eV. Chen 



10| et al. reported experimental and theoretical investigations of ZnS phonons in wurtzite 



and zincblende structure. They observe a good agreement between ab initio calculations 



and measurements. AGF is used in the literature to calculate ballistic phonon transport 
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- |l4| and results show good agreement with experiment. Other models that are used to 
calculate phonon transport are the acoustic mismatch model (AMM), which is only valid 
for the long wavelength phonons and the diffuse mismatch model (DMM), which assumes 
diffuse scattering at the interfaces 151]. Both methods do not take into account the structure 
of the interface and give therefore only an approximative description of phonon transport. 



II. METHOD 



The whole calculation divides up into two different tasks. The first task is the calculation 
of the equilibrium positions of the atoms and the determination of the interatomic force 
constants (ifcs). In the second task the AGF method is used to calculate the transmission 
function and the thermal conductance, whereas the ifcs from the first task are used to de- 
scribe the interatomic potential. With this approach only the harmonic part of the potential 
is considered, therefore anharmonic effects are not described. In principle anharmonic ef- 
fects can be incorporated in this approach by introducing additional self energies in the AGF 
method. 



A. Calculation of the Interatomic Force Constants 

In order to calculate the ifcs we use the code package ABINIT l| . This code is based on 
density functional theory and can be used with a pseudopotential method and a plane-wave 
expansion. For the exchange correlation potential, we use the local density approxima- 
tion(LDA). The pseudopotentials are represented in the TrouUier- Mart ins scheme 16|. The 
code uses a perturbation method to calculate the dynamical matrix on a discrete grid in the 
Brillouin zone. The use of a perturbation method avoids the use of super cells. Then the 
ifcs can be obtained by a discrete inverse Fourier transformation of the discrete dynamical 
matrix. Knowing the ifcs the dynamical matrix can be calculated in the whole Brillouin 
zone. Furthermore, also based on the same perturbation method the code is able to calcu- 
late the dielectric permittivity tensor and the born effective charges. Knowing the dielectric 
permittivity and the born effective charges the ifcs can be decomposed in a long range part 
which describes the dipole-dipole interaction and a short range part which describes the 
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electronic contribution to the interaction. For details see References 



17| and (iSj. 



B. Transport Calculation using Atomistic Green's Function Method 



A good overview of the AGF method applied to phonons is given in Reference 
Here we will summarize the basic features and show how to connect these method to the 
ifcs obtained from the ab initio calculation. First, we define the interaction of the atoms 
in the system and set up the harmonic matrix H. In the literature these interactions are 
often described by parameterized potentials {12!. 11^ . We use the ifcs obtained from ab inito 
calculations to describe the interaction between the atoms. The interaction between atom 
/i and atom 77 is represented by a 3x3 matrix. 

and are the masses of the atoms and k^^^j. is the force constant that describes 
the impact on atom /i in direction i if atom 77 is dislocated in direction j and vice versa. If 
we describe a structure that has periodicity in all 3 dimensions we can built the dynamical 
matrix which is the Fourier transform of the harmonic matrix. 

D{q) = Y,H{R)e-^'^^ (3) 

R 

i? is a lattice vector in real space and H{R) contains all interactions between cells that 
are connected by R. The phonon dispersion relation can be obtained by the square root of 
the eigenvalues of the dynamical matrix. The density of states can be calculated using the 
Green's function of the system. 

G{w,q) = iu'l-Diq))-' (4) 

The Green's function and the density of states are connected by an integration over the 
Brillouin zone 



n[LO) = 

JbZ 



f {Giu,q}-G{u,q}^)dq (5) 

JBZ 

For transport calculations we consider a system that is coupled to two semi-infinite con- 
tacts. The considered system is in general not periodic in transport direction. Therefore, we 



divide our system into layers of atoms that are perpendicular to the transport direction. This 
allows us to distinguish between interlayer and intralayer interactions. The layers are still 
periodic perpendicular to the transport direction. Consequently we can Fourier transform 
the interlayer and intralayer interactions along these directions. The periodicity in-plane is 
described by a 2-dimensional lattice with the lattice vectors i?^, where g is the layer index. 
The Fourier transform of the interlayer and intralayer interactions are given by 



H, 



E 



H{Rl)e-'^^^^ (intralayer) 



TgM) = Y,H{Ro)e~''^-'''' (interlayer) 

Ro 



(6) 



(7) 



Ro are the vectors that connect cells of neighbouring layers gg. With this 2-dimensional 
Fourier transformation we introduce a new Vector g^, a 2-dimensional unit cell, and the 
corresponding 2-dimensional Brillouin zone. 

The overall matrix describing the system has then the following form. 



T, 



9-^,9 
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(8) 



9+2,9+1 



Hg is the Fourier transformed intralayer interaction of layer g. Tg^g^i is the Fourier 
transformed interlayer interaction of layer g with g-|-l. In a typical calculation we take into 
account more than next nearest layer interaction. 

This matrix is infinite. The infinite matrix can be reduced to a finite matrix by using the 
concept of self-energies. The overall Green's function of the system is given by 



G{w, Qp) 



(uj^I -H -E, 



(9) 



where H is a finite submaxtrix of H(g^) and S/, and T,^ are the self energies for the left 
and right contact, respectively Uj. The self energies can be calculated using 
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Si = ngiTl (10) 

where Tl and Tn are the interlayer interactions which connect the middle region to the 
left and the right semi-infinite leads, gi and gr are the surface Green's functions of the left 
and right contact. We calculate these surface Green's functions using decimation techniques 
19|. 

From the self energies the broadening matrices can be obtained using 

Tl.K = t{^L,R - 4,if) (11) 

Consequently the transmission function of the system can be calculated using 



t{u,q,) = Tr[TLGTRG^] (12) 

The average transmission function per unit cell is given by an integral over the 2- 
dimensional Brillouin zone 



^^^^ ^ J2^ jj^z^^^'^^^^^^ ^'^'^^ 
The unit of this quantity is transmission per area. 

Knowing the average transmission one can calculate the total energy fiux per unit area 
J in the linear response 

J= — r hujt(uj)^^^^^dujAT (14) 

where f{uj,T) is the occupation function for the phonons. 
The conductance per unit area is defined as 



K 



J 1 ,dfico,T) 



hut{u) — — du (15) 



AT 2iv Jo ^ ' dT 

Bce con- 



If the exact structure of the interface is unknown one can estimate the inter 
ductance between material a and material b within the diffuse mismatch model 



15|. For 



the transmission function we assume that at the interfaces the phonons lose memory of 
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their original state and are either scattered in material a or material b. The corresponding 
probability is proportional to the transmission function of the bulk materials. The transmis- 
sion function is normalized such that for an interface between the same material the overall 
transmission is 1/2 of the bulk transmission. 

tiu) = ^(^^M^ (16) 
taioj) and ti,{uj) are the bulk transmission functions. 

We are using this model to estimate the influence of interfaces between two different 
materials. In future work we plan to investigate in addition the coherent transport across 
interfaces. 



III. RESULTS AND DISCUSSION 

Ground state calculations and cell shape optimization are performed using a k- 
point sampling of 6x6x3 for ZnO wurtzite and ZnS wurtzite and a cutoff energy of 60 hartree 
for the plane wave expansion. For the ZnS in zincblende structure we perform ground state 
and cell shape optimization with a k-point grid of 6x6x6 and a cutoff energy of 60 hartree 
for the plane wave expansion. Using these parameters the lattice constants are converged 
within 0.1%. The results of the cell shape optimization and results from other groups are 
shown in Table [H The lattice constants for ZnO wurtzite are about 1.6% smaller than the 
measured ones. The volume is about 3.9% too small. For ZnS in wurtzite structure and 
zincblende structure we also obtain lattice constant that are smaller than the measured 
one. This overestimation of the binding energies is a well known effect of the LDA. Based 
on the ground state density the dynamical matrix is calculated on a q-point grid in the 
brillouin zone. For wurtzite structure a 6x6x3 grid is used and for zincblende a 6x6x6 grid 
is used. The ifcs are calculated by a Fourier transform of the dynamical matrix. The ifcs 
are converged such that the conductance obtained from the ifcs are converged within 1%. 
The dispersion relation and density of states calculated with abinit are shown in Fig{T]for 
the three different material systems. The density of states and dispersion relation are in 



good agreement to calculations reported by other groups [10|, |20| . Figj2] left column shows 
the density of states of all three material systems. It shows that the overlap between the 
density of states of ZnO and ZnS is rather small regardless of the structure of ZnS, which 
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indicates a high phonon scattering at a ZnS/ZnO interface. We will emphasis this later in 
this paper by applying the diffuse mismatch model (DMM) to the two materials. 




FIG. 1: Dispersion relation and density of states. Top ZnO in wurtzite structure. Center ZnS in 
wurtzite structure. Bottom ZnS in zincblende structure 

Transmission functions for the pure materials in the different structures for different 
transport directions are shown in Figj2] middle column. The transmission functions of the 
wurtzite structures show a strong dependence on the transport direction. For ZnO the trans- 
mission function in c-direction is approximately a factor of 2-3 larger then for the other two 
directions. For ZnS the transmission function for the c-direction is also the highest. The 
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FIG. 2: Row 1: ZnO in wurtzite structure. Row 2: ZnS in wurtzite structure. Row 3: ZnS 
in zincblende structure. Column 1: Density of states. Column 2: Transmission function for 
different directions. Column 3: Temperature dependend conductance k. C, A and M stands for 
the directions pependicular to the respective c-,a-,m-plane. The Miller indices for these directions 
are [0001], [1120] and [OlIO] respectively 



transmission function for the zincblende structure shows no significant dependence on the 
direction. Using Eq.( IT5|) we calculate the conductance. In Figj3]the differential conductance 
which is in principle the integrand of Eq. f|T5|) is shown for different temperatures. It is shown 
that for low temperatures the main contribution of the conductance comes from phonons 
which have low energy. Raising the temperature leads to a small increase of the differential 
conductance at low energies and a strong increase at higher temperatures. At high temper- 
atures the shape of the differential conductance is only govern by the transmission function. 
Now we consider the conductance of the pure ideal materials. In particular, we consider 
only ballistic transport thus no scattering within the material. The only scattering is due to 
contact resistance which is also called Sharvin resistance. Therefore, the conductance shown 
in Figj2] has the unit of an interface conductance. The conductance for the pure materials 
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a(A) c(A) c/a u V(A3) 



ZnO Wurtzite 



This Work 3.198 



Theo. 



Exp. 



23 



Theo. 20] 3.198 



m 3.23 



Exp. 22 1 3.2499 



3.258 



5.167 1.615 0.379 22.882 

5.167 1.615 0.379 22.882 

5.168 1.6 0.381 23.3468 
5.20658 1.602 0.3819 23.8119 
5.220 1.6022 0.382 23.9924 



This Work 3.755 
Theo. 



Exp. 



IPC 



3.81 
24] 3.8227 



ZnS Wurtzite 

6.168 1.643 0.374 37.659 

1.33 1.64 - 39.2753 

6.2607 1.6378 0.3748 39.6154 



ZnS Zincblende 



This Work 5.32 
Theo. 



Exp. 



Id] 



25] 



5.4 
5.4109 



37.642 
39.366 
39.6049 



TABLE I: Calculated lattice constants compared with reported experimental and theoretical values 

in different transport directions are shown in Figj2] right column. ZnO has the highest con- 
ductance in C-direction. In the other two directions ZnO has a conductance that is reduced 
by a factor of 3 compared to the C-direction. This behavior is already included in the trans- 
mission function of ZnO. The conductances for ZnS shows only a small difference between 
the different directions in zincblende structure and the C-direction in wurtzite structure. A 
significant smaller conductance shows ZnS in wurtzite structure in A- and M-direction. The 
temperature dependence of all conductances show a similar behavior. At low temperatures 
the conductance increases rapidly with temperature until all phonon states are sufficiently 
occupied. Further increase in temperature leads only to a small increase in the conductance. 



The diffuse mismatch model is used to estimate the interface conductance of different 
combination of ZnO/ZnS interfaces. The conductance is again calculated using Eq. ffTSj) . The 
transmission function for an interface is calculated using Eq. (fT6|) . Since the transmission 
functions of ZnO and ZnS have no overlap above 28meV, all states above 28meV can not 
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FIG. 3: Differential conductance for ZnO in C-direction at temperatures of 100 K, 200K and 300K 

contribute to heat transport, which leads to a decrease in the conductance especially for 
higher temperatures. Since the transmission function for all directions in zincblende struc- 
ture for ZnS have a similar shape we used the 100 direction representative for all directions. 
FigJH shows the transmission function and conductance of the interface between ZnO and 
ZnS in 100 direction. Additionally we investigated interfaces between ZnO and ZnS where 
both are in wurtzite structure. Here we focus our investigations on interfaces with the same 
directions. Figj5] shows the transmission function and conductance of the interface between 
ZnO and ZnS both in wurtzite structure. These interface conductances are approximately 
one order of magnitude higher than typical interfaces conductances reported in the litera- 
ture 0|, which is reasonable, since the DMM gives an upper estimation for the interface 
conductance. The impact of such interfaces on the figure of merit is estimated by comput- 
ing the total conductance per m^ of a ZnO/ZnS superlattice structure with different period 
length. The total conductance per m^ is calculated using 



L is the length of the structure and n is the number of interfaces. KznO and KznS are the 
bulk conductivities and Kjnt is the interface conductance. Figl6] shows the figure of merit 
with different interfaces divided by the figure of merit without interfaces. In this calculation 
we assume that the interfaces have only an impact on the phonon thermal conductance 
and that the phonon thermal conductance is much larger than the electric part. Thus 
the ratio of {ZT)int with interfaces and by ZT without interfaces is given by nzr/i^zTint- 
For ajnt we used the average value of the interface conductance of ZnO/ZnS/wurtzite and 
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FIG. 4: Left: Transmission function of a ZnO/ZnS100(zincblende) interfaces using diffuse mis- 
match model (DMM) for different orientations of ZnO. Right: Corresponding thermal conductance. 
ZnOC-ZnSlOO stands for an interface between ZnO in C-direction and ZnS in 100 direction 
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FIG. 5: Left: Transmission function of a ZnO/ZnS(wurtzite) interfaces using DMM for different 
orientations of ZnO. Right: Corresponding thermal conductance 

ZnO/ZnS/zincblende presented in this paper. All parameters we used for this plot are listed 
in Table Ull It can be seen, that the figure of merit increasing linearly with the number 
of interfaces. The slope of this curve is given by the inverse of the interface conductance. 
Consequently superlattice structures of ZnO/ZnS are promising for a substantial increase of 
the figure of merit in thermoelectric devices. 



IV. SUMMARY AND CONCLUSION 



We performed DFT calculations of ZnO in wurtzite structure and ZnS in zincblende and 
wurtzite structure in the LDA. We calculated the lattice constants for this structures which 
agree well with other reported values. From DFT calculations we also obtain phonon band 
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FIG. 6: Left: Sketch of a ZnO/ZnS superlattice structure of length L. Right: Figure of merit as 
a function of the number of interfaces within the superlattice structure. The used parameter are 
given in Table HIl 



Parameter 



I^ZnO 
HZnS 



2S] 
28| 



I^ZnO/ZnSw 
I^ZnO/ZnSzb 

L 



Value 



54W/mK 
27W/mK 
0.175-109 W/m^K 
O.M-IO^ W/m^K 
1.0 • 10~'^m 

TABLE II: Parameters and references used for Fig. [6l The interface conductance nznO/ZnSw ^^i^d 
KznO/ZnSzb are the average values of the different interfaces shown in FigHJand FigjH The values 
are taken at a temperature of 300K. KznO is the ZnO bulk thermal conductivity at 300K average 



over the different directions 



28[ . KznS is the ZnS zincblende bulk thermal conductivity at 300K 



average over the different directions 



structures and phonon density of states which also agree well with reported ones by other 
groups. Based on these calculations we can also obtain the interatomic force constants (ifcs). 
These ifcs were then used within an AGF-Method to calculate bulk transmission functions. 
The difuse missmatch model (DMM) is used to estimate the interface conductance between 
ZnS and ZnO. We observed that the existence of such interfaces can have a significant impact 
on the thermal conductance and could lead to a substantial increase in the figure of merit. 
The DMM is only a rough method to estimate the interface conductance. In order to improve 
the description of the interface one can use the AGF method to built layer structures and 
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calculate the transmission function of an interface directly from Eq. (fT3|) . which we plan in 
the future. 



V. APPENDIX 

ZnO and ZnS are both intrinsic insulators which have non-vanishing effective charges, 
which means that the ifcs exhibits a dipol-dipol interaction. These dipol-dipol interaction 
shows an decrease in real space which makes it hard to converge the calculation by 
summing up all interactions in real space. To avoid this problem ABINIT decomposed 
the ifcs in a short range part and a dipol-dipol part. The dipol-dipol part is calculated in 



reciprocal space (for details see [18|). Unfortunately, in the AGF method the ifcs has to be 
defined in real space and therefore only a finite number of interactions can be taken into 
account. To estimate the error caused by this, we compared the dispersion relation and 
density of calculated with ABINIT and the dispersion relation and DOS calculate with the 
AGF. The DOS can be calculated in the AGF by using Eq.([n]). Fig|7] shows the dispersion 
relation calculated with the different methods. It can be seen, that both are almost equal 
except of three optical phonon modes which deviate in the nearness of the F-point. Fig|8] 
shows the DOS calculated with the different methods and it can be seen that there is only 
a very small difference due to the mentioned three optical phonon modes. Since there are 
only three optical phonon modes which are described wrong only around the F-point we 
can conclude that the error we made by truncating the interactions has only a small effect 
on total conductance. To estimate the radius after which we can truncate the sum of the 
interaction we perform convergence tests with respect to the transmission function. In FigH] 
the conductance for ZnO in C-direction is calculated using different radii for the interactions. 
After a radius of 1.12nm the conductance is convergend within 1%. The visible very small 
oscillations for large cutoff radii are due to the mentioned wrongly describes optical modes 
close to the F-Point. In this paper we use a cutoff radius of at least 1.1 2nm. 
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FIG. 7: Dispersion relation calculated with ABINIT (dashed) and the AGF method (solid) 




Energy [me V] 

FIG. 8: Desity of states calculated with ABINIT (dashed) and the AGF method (sohd) 




FIG. 9: Conductance of ZnO in C-direction calculated as a function cutoff radii at different tem- 
peratures 
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